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Following  our  integrated  hierarchical  modeling  framework  of  natural  gas  internal  reforming  solid  oxide 
fuel  cell  (IRSOFC),  this  paper  firstly  introduces  the  model  libraries  of  main  balancing  units,  including  some 
state-of-the-art  achievements  and  our  specific  work.  Based  on  gPROMS  programming  code,  flexible  con¬ 
figuration  and  modular  design  are  fully  realized  by  specifying  graphically  all  unit  models  in  each  level.  Via 
comparison  with  the  steady-state  experimental  data  of  Siemens-Westinghouse  demonstration  system, 
the  in-house  multi-level  SOFC-gas  turbine  (GT)  simulation  platform  is  validated  to  be  more  accurate  than 
the  advanced  power  system  analysis  tool  (APSAT).  Moreover,  some  units  of  the  demonstration  system  are 
designed  reversely  for  analysis  of  a  typically  part-load  transient  process.  The  framework  of  distributed 
and  dynamic  modeling  in  most  of  units  is  significant  for  the  development  of  control  strategies  in  the 
future. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

As  a  candidate  for  more  efficient  power  generation  than  the 
conventional  power  plants,  the  combined  cycle  by  solid  oxide  fuel 
cell  (SOFC)  and  gas  turbine  (GT)  attracts  extensive  attention  in  the 
worldwide.  Modeling  and  simulation  plays  an  important  role  in  the 
development  of  SOFC-GT  hybrid  generation  technology.  In  gen¬ 
eral,  various  system  configurations  and  numerical  requirements 
call  for  modular  simulation  with  different  complexities.  Recently, 
we  have  developed  an  in-house  multi-level  simulation  platform  for 
SOFC-GT  hybrid  generation  system  in  the  commercial  environment 
of  gPROMS  [1  ].  The  integrated  hierarchical  modeling  framework  of 
natural  gas  internal  reforming  SOFC  has  been  firstly  introduced  in 
our  previous  work  [2]. 

In  addition  to  the  multi-scale  modeling  of  fuel  cell  and  stack, 
models  of  other  balancing  units  including  the  gas  turbine,  reformer, 
heat  exchanger,  ejector,  burner,  splitter,  mixer,  etc.  are  indispens¬ 
able  for  system-level  analysis.  Table  1  exhibits  the  representative 
work  on  modeling  and  simulation  of  SOFC  or  molten  carbonate 
fuel  cell  (MCFC)  combined  heat  and  power  plant  and  SOFC-GT  or 
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MCFC-GT  hybrid  system  in  the  worldwide  institutes  [  1  -22  ].  In  gen¬ 
eral,  the  model-based  analyses  of  high  temperature  fuel  cell  power 
plant  can  be  summarized  as  follows: 


(1)  Analysis  of  full-load  and  part-load  operation  performance  and 
parameter  optimization  based  on  the  first  and  second  law  of 
thermodynamics. 

(2)  The  safe  or  healthy  operation  range  of  parameters  includ¬ 
ing  the  SOFC  temperature,  turbine  inlet  temperature  (TIT), 
pre-reforming  temperature,  surge  margin  of  the  compressor, 
steam-carbon  ratio,  and  pressure  difference  between  the  anode 
and  cathode  side  in  case  of  gas  recycling,  etc. 

(3)  The  necessity  of  introducing  controllers  and  various  control 
strategies,  for  example,  the  fixed  gas  flow  rate  or  fixed  gas 
utilization,  fixed  or  variable  turbine  speed,  the  different  tem¬ 
perature  control  modes  including  NFCRC  mode,  NETL  mode  and 
cascading  control,  etc. 

(4)  Technology  of  networked  stacks  or  stack  stage  which  focuses 
on  analysis  of  local  fuel  and  air  utilization,  current  and  tem¬ 
perature  distribution  among  the  serially  or  parallel  connected 
stacks.  And  some  new  configurations,  such  as  FC-Reformer-FC 
concept  by  using  an  external  reformer  as  the  intercooler  of  the 
upstream  stack  outlet  air. 
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Nomenclature 

A  Area  (m2)  or  pre-exponential  factor  in  Eq.  (12) 

a,  b,  c  Fitting  variables 
c  Concentration  ( mol  nrr 3 ) 

C  Heat  capacity  rate  (W  I<-1 ) 

cp  Mass  specific  heat  (J  kg-1 I<-1 ) 

d  Diameter  (m) 

E  Activation  energy  (J  mol-1 ) 

h  Mass  specific  enthalpy  (Jkg-1)  or  convective  heat 

transfer  coefficient  (J  m-2  s-1 ) 

Ah  Mass  specific  enthalpy  difference  or  polytropic  head 
(Jkg-1) 

k  Adiabatic  index 

I<e q  Equilibrium  constant  (bar2  or  1 ) 

M  Molecular  weight  (kg  mol-1 ) 

fn  Mass  flow  rate  (kg  s-1 ) 

N  Turbomachinery  rotational  speed  (rpm) 

Ncr  Corrected  speed  parameter  ( rpm  K-0-5 ) 

NTU  Number  of  heat  transfer  units 

p  Pressure  (Pa) 

p0  Standard  atmospheric  pressure  (Pa) 

P  Power  (kW) 

r  Volumetric  reaction  rate  (mol  m-3  s-1 ) or  radius  (m) 

R  Universal  gas  constant  (8.314 J  mol-1 I<-1 ) 

t  Time  (s) 

T  Temperature  (I<) 

u  Velocity  (ms-1) 

V  Volumetric  flow  rate  (m3  s-1 ) 

x  Molar  fraction 

z  Axial  coordinate  (m) 

Greek 

a  Total  heat  transfer  area  per  volume  of  heat 

exchanger  (m2  m-3) 

/3  parameter  in  (5  line  model 
s  Effectiveness  of  heat  exchanger 

£P  Void  of  packed  bed 

p  Efficiency 

k  Heat  conductivity  (W  m-1 1<-1 ) 

li  Dynamic  viscosity  (Pa  s) 

it  Pressure  ratio 

0  Corrected  mass  flow  parameter  (kg  s-1 1<0,5  bar-1 ) 

p  Density  (kg  m-3) 

cr  The  ratio  of  flow  area  (Afree)  to  frontal  area  (Afr)  of 

heat  exchanger  (m2  m-2) 

Subscripts  and  superscripts 
0  Design  point 

c  Compressor  or  cold  fluid 

cat  Catalyst 

g  Gas 

h  Hot  fluid 

in  Inlet 

max  Maximum 

min  Minimum 

out  Outlet 

ref  Reference  or  reforming  reaction 

s  Isentropic  process 

shift  Water  gas  shift  reaction 

t  Turbine  or  total 

w  Wall 


(5)  Methods  to  improve  the  flexibility  of  the  system  by  adopting 
an  external  combustor  and  bypass  valves  for  fast  start-up  and 
temperature  control. 

This  paper  firstly  introduces  the  model  library  of  main  balanc¬ 
ing  units  in  our  SOFC-GT  simulation  platform,  which  includes  some 
achievements  in  the  state-of-the-art  literatures  and  our  specific 
work.  Secondly,  in  a  gPROMS  commercial  environment,  all  unit 
models  in  different  level  are  graphically  specified  and  validated 
by  the  steady-state  experimental  data  of  Simense-Westinghouse 
demonstration  system.  Finally  a  transient  operation  with  variable 
load  is  analyzed  based  on  the  reverse  design  of  units. 

2.  Model  libraries  of  balancing  units 

2.1.  Turbomachinery 

2.1.1.  Isentropic  efficiency  model 

The  simplest  and  the  most  general  method  of  modeling  com¬ 
pressor  and  turbine  (C/T)  is  the  so-called  isentropic  efficiency 
model,  that  is 

_  fos,c(7s,c)  -  ^in(^in)  _  ^in(^in)  ~  hout(jout)  n 

/c  “  hout(Tout)  -  hin(Tin)  ’  /r  -  hin(Tin)  -  MrSit)  1  ; 

where  pc  and  pt  are  the  compressor  and  turbine  efficiency,  hin  and 
hout  are  the  mass  specific  enthalpy  of  inlet  and  outlet  gas,  hs  is  the 
mass  specific  enthalpy  for  isentropic  process,  Tin  and  Tou t  are  the 
inlet  and  outlet  gas  temperatures.  For  ideal  gas  in  an  isentropic 
process,  the  discharge  temperature,  Ts  is  related  to  the  pressure 
ratio  of  compressor  or  turbine  (ttc  =  Pout/Pin*  =  Pin/Pout) 

rs,c  =  rin7rc(fc“1)/k,  Ts,t=Tin7tt-^k-^k  (2) 

where  k  is  the  ratio  of  specific  heats  or  adiabatic  index. 

For  a  mass  flow  rate  m,  the  power  of  compressor  or  turbine  is 
the  enthalpy  difference  between  the  inlet  and  outlet  gas 

Pc  =  mc(hout  —  ^in)?  =  dit(h[n  —  hout )  (3) 

2.1.2.  Some  state-of-the-art  turbomachinery  models 

The  Process  Model  Library  (PML)  of  flow  transportation  in 
gPROMS  provides  a  fan  model  of  centrifugal  compressor  [23]. 
According  to  Fan  laws,  the  compressor  performance  between  the 
operating  point  and  the  design  point  (with  subscript  0)  is  related 
as 

V  _  N  Ah  _  /  N  \2 

^“iv  a ho~\N^)  (  ] 

where  V  is  the  volumetric  flow  rate,  N  is  the  compressor  rota¬ 
tional  speed,  Ah  is  the  mass  specific  enthalpy  difference  between 
the  discharge  and  suction  gas  (or  polytropic  head).  And  the  com¬ 
pressor  efficiency  is  considered  independent  on  the  compressor 
speed. 

The  Jensen  &  Kristensen  (J/I<)  model  has  been  widely  used  in 
modeling  of  the  turbocharger  of  internal  combustion  engines  [24]. 
In  J/K  compressor  model,  the  dimensionless  head  parameter  and 
compressor  efficiency  are  expressed  as  functions  of  the  normalized 
flow  rate  and  the  inlet  Mach  number.  In  J/K  turbine  model,  in  addi¬ 
tional  to  the  nozzle  equation,  the  turbine  efficiency  is  expressed 
as  a  polynomial  function  of  the  corrected  turbine  speed  parameter 
and  the  blade-speed  ratio,  which  is  also  suitable  for  variable  nozzle 
turbines. 

Without  physical  meaning  of  ft,  line  model  is  a  purely  math¬ 
ematical  method.  As  shown  in  Fig.  1,  lines  are  a  set  of  equally 
spaced  parabolas  or  straight  lines,  which  have  a  unique  intersec¬ 
tion  point  with  each  constant  speed  line  [12,17,25].  By  setting  two 
parabolic  lines  at  /3  =  0  and  fi  =  \  (each  determined  by  three  data 
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Table  1 

Representative  work  on  modeling  and  simulation  of  SOFC/MCFC  based  system. 


Institute 

Country 

Main  research  contents  and  achievements 

National  Fuel  Cell  Research  Center  (NFCRC) 

U.S.A. 

Advanced  power  systems  analyses  tools  (APSAT),  energy  and  exergy  analyses  of  single-stage 
and  two-stage  SOFC-Humid  Air  Turbine  (HAT)  cycle,  necessity  of  introducing  controllers, 
different  temperature  control  mode  including  NETL  mode,  NFCRC  mode  and  cascading  control 
[3,4] 

Reasonable  air  recycling  in  the  cathode  and  technology  of  stack  stage,  and  comparison  among 
the  thermodynamic  performance  of  various  external  reforming  and  internal  reforming  hybrid 
configurations  [5] 

National  Energy  Technology  Laboratory 
(NETL) 

U.S.A. 

University  of  Wisconsin-Madison 

U.S.A. 

EES-based  modeling  of  SOFC  system  for  small-scale  stationary  application,  configuration  with 
anode  and  cathode  gas  recycling,  optimal  system  design  and  comparison  among  three  control 
modes  [6] 

The  Pennsylvania  State  University,  Park 
University 

U.S.A. 

Lumped  and  control-oriented  modeling  of  MCFC  power  plant,  decentralized  controller  design 
for  the  stack  temperature,  pressure  difference,  fuel  utilization  and  steam-carbon  ratio  [7] 

University  of  Genova 

Italy 

TRANSEO  simulation  tool,  full-load  and  part-load  performance  of  SOFC/MCFC-GT  system, 
effects  and  control  strategies  of  gas  flow  rates,  gas  recirculation  ratio,  turbine  speed,  bypass 
valve,  etc.  [8,9] 

University  of  Perugia 

Italy 

Temperature  control  strategies  using  four-way  valve  in  an  atmospheric  MCFC-GT  power  plant 
[10] 

Modeling  of  combined  SOFC  CHP  plant  based  on  Pro/II,  analysis  of  various  parameters 
including  internal  reforming,  air  inlet  temperature,  fuel  utilization,  and  gas  recycling  in  the 
anode  and  the  cathode  [11] 

Juelich  Research  Center 

Germany 

Norwegian  University  of  Science  and 
Technology 

Norway 

gPROMS-based  toolbox  for  SOFC-GT  power  plant,  system  design,  steady-state  part-load 
performance,  transient  operation  and  decentralized  controller  [12] 

Lund  University 

Sweden 

AspenPlus-based  analysis  of  full-load  operation  of  SOFC-GT  system,  networked  SOFC  stacks 
[13,14] 

Delft  University  of  Technology 

Netherlands 

MCFC  power  plant  with  gas  recirculation  and  waste  heat  utilization  in  the  form  of  an  expander 
and  steam  production  cycle,  influence  of  operating  temperature  on  system  performance  [15] 

University  of  Waterloo 

Canada 

AspenPlus-based  modeling  and  simulation  of  Siemens  Westinghouse  SOFC-GT  hybrid  system 
[16] 

Component  matching  in  IRSOFC-GT  power  plant,  performance  analysis  of  heat  and  power 
cogeneration,  configuration  with  an  afterburner  between  SOFC  stack  and  GT  [17] 

Nanyang  Technological  University 

Singapore 

University  of  Tokyo,  Nippon  Telegraph  and 
Telephone  Corp. 

Japan 

Conceptual  design  of  a  30  kW-scale  SOFC-MGT  hybrid  system,  stack  stage  [18,19] 

Seoul  National  University 

Korea 

Various  part-load  control  modes  for  SOFC-only  system  and  SOFC-GT  hybrid  system  [20] 

Tsinghua  University,  University  of  Science 
and  Technology  Beijing 

PR  China 

AspenPlus-based  modeling  and  thermodynamics  analysis  of  SOFC-MGT  hybrid  system, 
multi-level  simulation  platform  of  SOFC-GT  hybrid  system  based  on  gPROMS  environment 
[1,2,21] 

Shanghai  Jiao  Tong  University 

PR  China 

Dynamic  modeling  and  part-load  performance  of  atmospheric  MCFC-GT  hybrid  system  [22] 

points),  all  the  polynomial  coefficients  of  ft  lines  can  be  determined. 
The  ft  lines  are  usually  selected  along  the  direction  of  the  high¬ 
est  efficiency  curve  (the  polygonal  line  in  Fig.  1).  Thus  the  flow 
rate,  pressure  ratio  and  efficiency  can  be  expressed  as  nonlinear 
polynomials  of  the  rational  speed  and  ft 


n1  n2  ”3 

<t>  =  J2  ayN^,  Jt=  byN^,  r]  =  ci]K tP  (5) 

1=0, j=0  i=0j=0  i=0,j=0 

where  0  =  m  (Tin)1/2/(pin/p0)  is  the  corrected  mass  flow  rate 
parameter,  Ncr  =  N/(Tin)1/2  is  the  corrected  speed  parameter,  Ncr  = 


NCr/Ncro  is  the  dimensionless  corrected  speed  parameter,  ni-n3 
are  orders  of  the  polynomials,  ay,  by  and  Cy  are  fitting  parameters. 

The  turbine  can  be  model  similarly  using  ft  line  grids.  However, 
the  turbine  inlet  chokes  above  a  certain  pressure  ratio,  i.e.  the  mass 
flow  is  no  longer  a  function  of  the  pressure  ratio  respectively  the 
parameter  ft.  The  following  elliptic  approach  instead  of  a  polyno¬ 
mial  function  is  preferred  to  describe  this  behavior  of  constant  mass 
flow  rate  [12] 

1/2 

,  z  =  c0  +  c1Ncr  +  c2JV;?r  (6) 

where  a,  b  and  c0-c2  are  fitting  parameters. 


cj)  =  a 


1  - 


2.1.3.  Neural  network  model 

Like  the  above-mentioned  J/I<  model  and  /3  line  model,  tur¬ 
bomachinery  modeling  is  usually  based  on  performance  map  by 
fitting  the  nonlinear  relationships  between  Ncr,  0,  tt,  and  p.  In  gen¬ 
eral,  two  types  of  functions  $=f{7t,Ncr)  and  tt  =/(</>,  Ncr)  are  both 
required  for  system-level  analysis,  especially  for  dynamic  simula¬ 
tion.  From  typical  speed  contours  of  compressor  as  shown  in  Fig.  1, 
close  to  the  surge  zone,  the  constant  speed  line  is  almost  parallel 
to  the  axis  of  flow  rate.  When  function  in  the  form  0=/(7r,Ncr)  is 
taken,  a  small  disturbance  of  pressure  ratio  leads  to  a  large  varia¬ 
tion  of  flow  rate,  even  a  multi-value  problem  appears.  Close  to  the 
choke  zone,  function  in  the  form  tt  =/(0,  Ncr)  encounters  the  similar 
problem. 

Although  p  line  model  solves  this  problem  in  some  extent,  it 
generally  requires  all  the  operating  data  should  be  within  the  enve¬ 
lope  curve  enclosed  by  the  parabolic  lines  /3  =  0  and  /3  =  1,  which 
means  a  poor  ability  of  extrapolation.  Neural  network  with  a  typical 


Fig.  1.  Typical  speed  contours  of  compressor  with  j3  line  grids. 
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Fig.  2.  Flow  chart  of  neural  network  model  of  turbomachinery.  Whereas,  r  is  the  geometric  ratio  of  the  actual  machine  to  the  original  machine  (superscript  as ')  for  similarity, 
and  6  is  the  rotation  angle  for  coordinate  transformation. 


topology  of  feed-forward  multilayer  perceptron  including  biases,  a 
sigmoid  layer,  and  a  linear  output  layer  is  capable  of  approximating 
any  function  with  a  finite  number  of  discontinuities. 

To  alleviate  the  problem  of  sharp  slope  in  C/T  map,  the  tech¬ 
nology  of  coordinate  transform  or  Givens  transform  is  firstly 
introduced  for  map  preprocessing  [26].  During  network  train¬ 
ing,  some  preprocessing  technologies,  such  as  data  normalization 
and  early  stopping  are  used  to  improve  the  network  general¬ 
ization.  For  a  general  design  based  on  a  same  original  map,  the 
similarity  theory  is  also  utilized  for  the  size  match  [26].  Via  train¬ 
ing  the  weights  and  biases  of  network,  a  good  interpolation  and 
extrapolation  accuracy  can  be  obtained  from  network  output  after 
postprocessing  of  inverse  coordinate  transform  and  inverse  simi¬ 
larity.  Fig.  2  shows  the  flow  chart  of  neural  network  model  in  our 
work.  For  the  surge  line  of  centrifugal  compressor  map,  the  pressure 
ratio  can  be  simply  expressed  as  a  three-order  polynomial  of  the 
flow  rate. 

Considering  the  effective  inertia  of  shaft,  all  the  above  C/T  mod¬ 
els  (except  for  isentropic  efficiency  model)  can  be  easily  expanded 
for  dynamic  analysis. 

2.2.  Steam  methane  reformer 

Similar  to  the  multi-level  SOFC  modeling  [2],  the  model 
library  of  steam  methane  reformer  also  includes  quasi-equilibrium, 
lumped  dynamic  and  distributed  dynamic  models.  Consid¬ 
ering  the  steam  reforming  reaction  (CFI4  +  H20^>C0  +  3Fl2), 
water-gas  shift  reaction  (C0  +  H20^>C02 +  H2)  and  overall  reac¬ 
tion  (CFI4  +  FI20^>C02  +4FI2)  in  a  fixed  bed  with  Ni-based  catalyst, 
thermodynamic  equilibriums  are  assumed  in  the  quasi-equilibrium 
model 


^q,ref*CH4*H20  =  *H2xCO,  ^eq,shift*CO*H20  =  *H2*C02  (7) 

where  x*  is  the  molar  fraction  of  species  i,  p  is  the  total  pressure, 
I<e q  is  the  equilibrium  constant. 

The  standard  Langmuir-FIinshelwood  expression  of  reaction 
kinetics  is  introduced  in  lumped  and  distributed  dynamic  modeling 
[27] 


kref 

Fref  =  A:at  ^5 
PH2 


PCH4Ph20  (Pj-^PCo/^eq.ref) 

(DEN)2 


rshift  —  Pc  at 


fcshift  PcoPh2o 
Ph2 


(Ph2PC02  /^eq.shift ) 
(DEN)2 


2 

_ ^  ^overall  ^H4Ph2o 

'overall  Pc  at  ^  5 

Ph2 


(Pll2Pc  02  /fQq,  re  f^eq,  shift) 

(DEN)2 


DEN  =  1  +  KcqPco  +  Kh2Ph2  +  Kch4Pch4  +  KHfPH2° 


(8) 

(9) 

(10) 

(11) 


where  pc at  is  the  catalyst  density,  rref,  rshift  and  roverall  are  the 
reaction  rates  of  reforming,  water  shift  and  overall  reaction,  respec¬ 
tively. 

Table  2  lists  the  kinetic  constants  X=kref,  /<Shift»  ^overall  or  JQ  in 
Eqs.  (8)— (1 1 )  in  Arrhenius-type  expression  [27] 

x=Aexp(-T)(T__T)  (12) 

For  distributed  modeling  of  a  tubular  packed  bed  reactor,  pres¬ 
sure  drop  along  the  axial  coordinate  can  be  described  by  Ergun 
equation  [23],  and  the  convective  heat  transfer  coefficient  between 
gas  and  solid  phase,  hg  can  be  calculated  empirically  by  [28] 


where  dcat  is  the  catalyst  diameter,  rin  is  the  inner  radius  of  tube,  £p 
is  the  void  of  packed  bed,  pg,  /xg,  /cg  and  ug  are  the  density,  dynamic 
viscosity,  heat  conductivity  and  velocity  of  gas. 

By  neglecting  the  dynamics  of  total  gas  pressure,  3p/3t  =  0  and 
3ct/3t=  d{plRT)ldt=  -(ct/T)  3l/3t,  the  total  gas  mass  balance  is  gov¬ 
erned  by  an  ordinary  differential  equation  instead  of  a  partial 
differential  equation  [2].  Although  this  treatment  brings  stronger 
numeric  couple  between  mass  balance  and  energy  balance  of  gas 
phase,  it  is  more  reasonable  than  the  assumption  of  constant  gas 
velocity  in  some  literatures  [12,23]. 

2.3.  Heat  exchanger 

The  £-NTU  method  is  the  most  general  for  thermodynamic  mod¬ 
eling  of  heat  exchangers.  Some  analytical  relationships  between  the 
effectiveness  (£)  and  number  of  heat  transfer  units  (NTU)  of  heat 
exchanger  in  different  flow  modes  [29]  are  collected  in  our  work. 

For  counter  flow  which  is  most  commonly  used  in  tube-shell 
heat  exchangers,  the  £-NTU  function  is 

1  -  exp[-NTU(l  -  (Cmin/Cmax))] 

1  -  (Cmin/Cmax)  exp[-NTU(l  -  (Cmin/Cmax))]  1  j 

where  Cmin  and  Cmax  means  the  minimum  or  maximum  heat  capac¬ 
ity  rate  of  hot  or  cold  fluids.  When  Cmin  =  Cmax,  £  =  NTU/(1  +  NTU). 


Table  2 

Kinetic  parameters  in  Eqs.  (8)— (1 1 )  by  Arrhenius-type  expression. 


X 

A 

E  (kj  mol-1) 

Tref(K) 

kref  (kmol  bar0-5  kg^r1  hr1 ) 

1.842  x  10-4 

240.1 

648 

kshift  (kmol  bar-1  kg^r1  hr1) 

7.558 

67.13 

648 

^overall  (kmol  bar0-5  kgcat_1  h-1) 

2.193  xlO-5 

243.9 

648 

Kh2  (bar-1) 

0.0296 

82.9 

648 

Xh2o 

0.4152 

88.68 

823 

KCh4  (bar-1) 

0.179 

38.28 

823 

I<co  (bar-1) 

40.91 

70.65 

648 
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For  cross  flow  of  two  unmixed  fluids,  which  is  often  configured 
in  some  compact  heat  exchangers,  such  as  the  recuperator  of  gas 
turbines  and  vehicular  radiators,  it  is 


s  =  1  -  exp 


exp  [-NTU(Cmin/Cmax)NTU-0-22]  -  1 


NTU 


r-0.22 


(15) 


In  our  previous  work,  we  have  presented  the  quasi-static  rating 
and  sizing  method  of  compact  heat  exchangers  in  fuel  cell  systems 
[30].  We  further  develop  a  distributed  dynamic  model  with  a  set 
of  uniform  governing  equations  for  both  tube-shell  and  compact 
heat  exchangers. 

As  for  concurrent  or  countercurrent  flow,  along  the  axial  coor¬ 
dinate  z,  the  energy  balance  of  cold  (k  =  c)  or  hot  (k  =  h)  fluid  is 

Pkcp,k^-  =  -Pkukcp,k-fif  +  ~l/o,khk(Iw  -  7k)  (k  =  C,  h)  (16) 

where  p ,  cp,T,u  are  the  fluid  density,  mass  specific  heat,  temper¬ 
ature  and  velocity,  Tw  is  the  wall  temperature,  h  is  the  convective 
heat  transfer  coefficient,  a  is  the  total  heat  transfer  area  per  volume 
of  heat  exchanger,  o  is  the  ratio  of  free  flow  area  (Afree)  to  frontal 
area  (Afr)  at  given  fluid  side,  rj0  is  the  overall  surface  temperature 
efficiency. 

For  gaseous  fluid,  using  the  ideal  gas  law  ( p=pM/RT )  and 
neglecting  the  pressure  dynamics  (3p/3t  =  0),  the  mass  balance  of 
given  fluid  is  governed  by 


0  =  -9(Pk"lJk)  +  ^%0.khkCTw  -  Tk) 
oz  Cp,k  &k 


(17) 


The  first  item  can  be  further  reduced  to  3(puT)/3z  =  ( pM/R)3u/3z 
under  assumption  of  constant  pressure.  For  incompressible  fluid, 
the  above  equation  is  replaced  by  a  uniform  velocity  distribution, 
3u/3z  =  0. 

The  energy  balance  of  wall  is  governed  by 


dTw  _  ^  3 2Tw  (  1  a  h  rT  t  \ 

PwCp.w  —^r  —  KW  ^2  ^  /  ^^free,k~rJo,knklIk  ~  *w) 


(18) 


k=c,h 


where  kw  is  the  heat  conductivity  of  wall,  Aw  is  the  section  area  of 
wall.  In  general,  the  adiabatic  boundaries  are  assumed. 

Without  considering  effects  of  fins, 


3o,k  —  1  ? 


Ok 

°k 


1 

rhyd,k 


(19) 


where  rhyd  is  the  hydraulic  radius  at  given  fluid  side,  Eqs.  (16)-(18) 
will  reduce  to  the  governing  equations  for  tube-shell  heat  exchang¬ 
ers. 


2.4.  Ejector 

For  a  gas  ejector,  the  key  geometries  are  the  minimum  section 
area  of  nozzle  and  the  section  area  of  mixing  chamber,  and  the 
representative  thermodynamic  parameters  include  the  mass  flow 
rate,  pressure  and  temperature  of  working  flow,  induced  flow  and 
compressed  flow  [31].  Governed  by  mass,  momentum  and  energy 
balance,  the  quasi-static  gas  ejector  model  for  anode  recirculation 
in  fuel  cell  systems  has  been  presented  in  our  previous  work  [32]. 
In  addition,  three  critical  conditions  are  also  considered  here,  i.e. 
the  induced  flow  velocity  reaches  the  critical  adiabatic  velocity  (or 
local  sound  velocity)  at  the  inlet  or  an  intermediate  section  of  mix¬ 
ing  chamber,  or  the  compressed  flow  velocity  reaches  the  local 
sound  velocity  at  the  mixing  chamber  outlet.  For  gas  recirculation 
in  fuel  cell  systems,  the  gas  ejector  model  can  be  further  simplified 
by  considering  the  induced  flow  and  compressed  flow  as  incom¬ 
pressible  fluids  due  to  small  pressure  difference  between  them 
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Fig.  3.  Graphic  specification  of  distributed  SOFC  model  on  the  tag  of ‘surface  diffu¬ 
sion’. 


2.5.  Others 

The  thermodynamic  models  or  continuous  stirred-tank  reactor 
(CSTR)  models  of  other  balancing  units,  such  as  burner,  splitter, 
mixer,  shaft,  pipeline  and  power  load,  are  also  included  in  our 
SOFC-GT  simulation  platform.  All  the  unit  models  in  each  level 
are  configured  with  a  graphic  specification  in  gPROMS  environ¬ 
ment.  Fig.  3  shows  an  example  of  configuration  for  distributed  SOFC 
model  on  the  Tag  of ‘surface  diffusion’.  In  additional  to  the  built-in 
stream  interfaces  (PMLMaterial  and  PMLControl)  in  gPROMS,  some 
special  interfaces  (radiation  interface,  power  load  interface)  are 
also  defined  for  high-level  configurations.  Thus,  users  can  deploy 
their  simulation  flexibly  for  modular  design  and  analysis. 


3.  Model  validation 

The  multi-level  unit  model  libraries  are  validated  here  by 
the  experimental  performance  of  Siemens-Westinghouse  Power 
Corporation  (SWPC)  220  kW  pressurized  tubular  SOFC-GT  demon¬ 
stration  hybrid  system  [33].  The  desulfurized  fuel  (natural  gas)  is 
fed  into  the  stack  to  be  reformed  and  react  electrochemically  with 
the  preheated  air.  The  cathode  outlet  gas  and  the  unrecycled  anode 
outlet  gas  continue  to  burn  in  the  post-combustion  zone  of  TSOFC. 
The  combustion  exhaust  firstly  preheats  the  compressed  air  via  the 
air  supply  tube  (AST),  then  goes  out  of  stacks  and  drives  a  high- 
pressure  (TIP)  turbine  to  compress  air  and  sequentially  drives  a 
low-pressure  (LP)  turbine  to  generate  power.  The  internal  reform¬ 
ing  inside  SOFC  stacks  includes  an  indirect  internal  reforming  (HR) 
process  in  an  external  reformer  and  a  direct  internal  reforming 
(DIR)  process  in  the  SOFC  anode.  Anode  recirculation  via  two  ejec¬ 
tors  is  used  to  provide  steam  in  the  external  reformer,  where  the 
overall  strongly  endothermic  reaction  is  supported  by  the  radia¬ 
tion  energy  from  SOFC  stacks.  Some  bypass  valves  are  deployed  for 
temperature  modulations.  Fig.  4  shows  the  gPROMS  realization  of 
the  system  configuration. 

Table  3  lists  the  characteristics  of  unit  models  in  Fig.  4.  For  a 
simple  simulation  of  full-load  and  stationery  operation,  most  of 
unit  models  are  quasi-static  in  thermodynamics.  Flowever,  a  CSTR 
model  for  the  external  reformer  and  a  distributed  dynamic  SOFC 
model  with  the  detailed  radiation  model  are  used  for  better  pre¬ 
diction  [2,34].  Because  of  lacking  of  the  C/T  performance  map,  the 
isentropic  efficiency  model  is  applied  in  this  section  for  all  the  tur¬ 
bomachineries,  and  £-NTU  model  is  used  for  the  recuperator  due  to 
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Pipe003 


Fig.  4.  gPROMS  realization  of  the  configuration  of  SWPC  SOFC-GT  demonstration  system. 


the  same  reason.  Pipelinel  -3  are  configured  to  describe  the  system 
heat  loss. 

Table  4  shows  the  structural  parameters  and  physical  prop¬ 
erties  of  SOFC.  The  interface-type  positive-electrolyte-negative 
(PEN)  transport  model  is  used  by  considering  the  small  anode 
thickness,  the  ohmic  polarization  is  calculated  by  transmission-line 
model,  and  the  activation  polarization  is  obtained  by  Arrhenius- 
type  expressions  of  activation  resistance.  Table  5  shows  some 
parameters  related  to  electrochemical  kinetics  [35]  and  Table  6 
shows  the  parameters  related  to  flow  and  heat  transfer  [12,36]. 

Tables  7  and  8  compare  the  system  performance  and  status 
parameters  of  stream  temperature  (°C),  pressure  (bar)  and  mass 
flow  rate  (kgs-1)  in  the  base  case  among  the  experimental  data, 
simulation  results  of  advanced  power  system  analysis  tool  (APSAT) 


Table  3 

Name  and  property  of  unit  models  in  Fig.  4. 


Unit  model  name  Model  function  and  property 


Fuel 

Air 

SinkOOl 

Stack 

Compressor 

HPTurbine 

LPTurbine 

Reformer 

Recuperator 

HTX001 

InjectorOOl 

Burner 

Load 

Splitter001-003 

Mixer001-002 

Shaft 

Pipe001-006 
Pipelinel -3 


‘Source’  model  in  gPROMS  for  inlet  natural  gas 

‘Source’  model  in  gPROMS  for  inlet  air 

‘Sink’  model  in  gPROMS  for  atmosphere 

Distributed  dynamic  model  for  tubular  SOFC  stack 

Isentropic  efficiency  model  for  compressor 

Isentropic  efficiency  model  for  high-pressure  turbine 

Isentropic  efficiency  model  for  low-pressure  turbine 

Lumped  dynamic  model  for  DIR  reformer 

e-NTU  thermodynamic  model  for  recuperator 

£-NTU  model  for  air  supply  tube  (AST)  in  TSOFC 

Gas  ejector  model  for  two  anode  ejectors 

Quasi-static  model  for  post-combustion  zone  in  TSOFC 

Power  load  model  for  generator 

Splitter  model  for  bypass 

Mixer  model  for  stream  mixer 

Power  balance  model  for  C/T  shaft 

Linear  flow  resistance  model 

Linear  flow  resistance  model  with  heat  loss 


Table  4 

Structural  parameters  and  physical  properties  of  TSOFC. 


Description 

Value 

Anode  thickness  (m) 

1.0  x  10-4 

Cathode  thickness  (m) 

2.2  x  10“3 

Electrolyte  layer  thickness  (m) 

4.0  xl0“5 

Anode  porosity 

0.4 

Anode  tortuosity 

3 

Cathode  porosity 

0.5 

Cathode  tortuosity 

1.5 

Electrolyte  conductivity  (S  irr 1 ) 

3.34  x  104  exp(-10300/T) 

Anode  electronic  conductivity  (S  irr 1 ) 

9.5  x  107/T-exp(— 1150/7) 

Cathode  electronic  conductivity  (S  nrr 1 ) 

4.2  x  107/T  exp(-1200/T) 

Interconnector  conductivity  (S  m-1 ) 

9.3  x  106/Texp(-1100/T) 

Cell  tube  outer  radius  (m) 

1.1  x  10“2 

Cell  length  (m) 

1.5 

AST  inner  radius  (m) 

2.5  x  10-3 

AST  outer  radius  (m) 

4  x  10-3 

Gap  between  cell  tubes  (m) 

2  x  10-3 

Interconnector  perimeter  (m) 

6xl0“3 

Interconnector  thickness  (m) 

1  x  10“3 

Serial  cell  number 

384 

Parallel  cell  number 

3 

Anode  density  (kg  m-3 ) 

7740 

Anode  mass  specific  heat  (J  kg-1  K-1 ) 

300 

Cathode  density  (kg  ur3 ) 

7740 

Cathode  mass  specific  heat  (J  kg-1 1<-1 ) 

300 

Anode  density  (kg  m-3 ) 

5300 

Anode  mass  specific  heat  (J  kg-1  K-1 ) 

400 

Electrolyte  density  (kg  m-3 ) 

6000 

Electrolyte  mass  specific  heat  (J  kg-1  K_1 ) 

500 

AST  density  (kg  m-3 ) 

3970 

AST  mass  specific  heat  (J  kg-1 1<-1 ) 

1167 

Anode  thermal  conductivity  (W  m- 1  K-1 ) 

6.23 

Cathode  thermal  conductivity  (W  m- 1  K-1 ) 

9.6 

Electrolyte  thermal  conductivity  (W  m-1 1<-1) 

2.7 

Electrolyte  thermal  conductivity  (W  m_  1  K-1 ) 

2.7 

PEN  emissivity 

0.8 

AST  emissivity 

0.2 
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(b)  Dimensionless  cell  length 


(c)  Dimensionless  cell  length 


Fig.  5.  Axial  distribution  of  (a)  temperatures  of  anode  gas,  cathode  gas,  preheated  air,  PEN  and  AST,  (b)  fuel  concentrations,  and  (c)  polarizations  in  SOFC. 


[33]  and  our  model.  In  our  calculation,  the  turbine  generation 
power,  SOFC  cell  voltage  and  SOFC  operating  pressure  are  set  as  the 
experimental  data.  It  is  considered  more  reasonable  to  get  the  aver¬ 
age  SOFC  temperature  from  the  temperature  distribution  instead  of 
setting  it  as  APSAT  did.  Here,  the  stack  operating  temperature  is  the 
average  value  among  axial  distribution  of  PEN  temperature,  and  the 
system  efficiency  is  calculated  based  on  the  lower  heating  value  of 
natural  gas  (49,81  OkJ kg-1 )  [33].  Mainly  due  to  our  advanced  dis¬ 
tributed  SOFC  model  with  the  detailed  radiant  model  [2],  most  of 
our  calculation  obtained  a  smaller  relative  error  (less  than  1%)  than 
what  APSAT  did. 

More  importantly,  superior  to  APSAT,  the  distributed  SOFC 
model  provides  the  information  of  local  parameters.  Fig.  5(a)  shows 
the  temperature  distributions  of  the  anode  gas,  cathode  gas,  pre¬ 
heated  air,  PEN  and  AST  along  the  tube  length.  Because  of  the 
strong  endothermic  reaction  of  methane  reforming,  both  the  gas 
temperature  and  solid  temperature  decrease  firstly  close  to  the 


Table  5 

Parameters  related  to  electrochemical  kinetics  in  TSOFC. 


Parameter 

Symbol 

Value 

Anode  activation  energy  (J  mol-1  )a 

Ea 

1.1  X 

105 

Cathode  activation  energy  (J  mol-1  )a 

Ec 

1.6  x 

105 

Pre-exponent  factor  of  H2  kinetics  (Am-2)3 

1<h2 

2.13 

x  108 

Pre-exponent  factor  of  CO  kinetics  (Am-2) 

1<c  0 

5.0  x 

107 

Pre-exponent  factor  of  02  kinetics  (Am-2)a 

ko2 

1.49 

x  1010 

H2  and  CO  reaction  order3 

m 

0.25 

02  reaction  order3 

m 

0.25 

anode  inlet,  then  increase  gradually  due  to  the  exothermic  elec¬ 
trochemical  reaction.  The  gap  between  the  highest  and  lowest  cell 
temperature  (almost  130 1<)  means  a  significant  temperature  gra¬ 
dient,  i.e.  the  temperature  of  hot  spots  is  much  higher  than  the 
average  temperature,  which  supports  the  necessity  of  distributed 
SOFC  modeling  with  reasonable  numeric  cost.  On  the  other  hand, 
at  the  anode  and  cathode  outlet  where  the  cell  temperature  is  the 
highest  under  co-flow  condition,  the  gas  temperatures  are  almost 
equal  to  the  local  solid  temperature  due  to  full  convective  heat 


Table  6 

Parameters  related  to  flow  and  heat  transfer  in  TSOFC. 


Parameter 

Symbol 

Value 

Nusselt  number  for  anode  gas  channel3 

Nua 

Zhukauskas 

correlation 

Nusselt  number  for  cathode  gas  channel3 

Nuc 

4.4 

Nusselt  number  for  inner  AST3 

NUAsT,in 

Gnilinski 

correlation 

Nusselt  number  for  outer  AST3 

Nu  AST, out 

6.0 

Friction  factor  of  anode  gas3 

/a 

80 

Friction  factor  of  cathode  gas3 

fc 

90 

Friction  factor  of  air  in  AST3 

/air 

64 

Pressure  loss  coefficient  at  anode  inlet5 

£a.in 

0.5 

Pressure  loss  coefficient  at  anode  outlet5 

£a,out 

0.5 

Pressure  loss  coefficient  at  cathode  inlet5 

£c,in 

0.5 

Pressure  loss  coefficient  at  cathode  outlet5 

£c,out 

0.5 

Pressure  loss  coefficient  at  AST  inlet5 

^AST.in 

0 

Pressure  loss  coefficient  at  AST  outlet5 

^AST.out 

1.0 

Radiation  heat  from  SOFC  to  pre-reformer  (kW) 

Q 

20 

aFrom  Ref.  [36]. 
bFrom  Ref.  [12]. 


From  Ref.  [35]. 
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Table  7 

Comparison  of  system  performance  among  experiments,  APSAT  and  our  model. 


Performance 

Experiments 

APSAT 

Our  model 

Turbine  generation  power  (kW) 

21 

21.6 

21 

SOFC  operating  temperature  (°C) 

995 

1000 

998.05 

SOFC  operating  pressure  (bar) 

2.9 

2.94 

2.9 

SOFC  cell  voltage  (V) 

0.639 

0.633 

0.639 

SOFC  stack  voltage  (V) 

244 

243.1 

245.376 

SOFC  stack  current  (A) 

700 

694.2 

699.65 

SOFC  stack  power  (kW) 

170.8 

168.78 

171.677 

Total  AC  power  (kW) 

183.45 

181.94 

184.09 

System  efficiency  (%) 

52.44 

52.0 

52.62 

transfer.  As  a  result,  the  gas  outlet  temperature  of  SOFC  is  actually 
much  higher  than  that  from  lumped  modeling,  which  is  meaningful 
for  the  control  strategy  of  TIT. 

Fig.  5(b)  shows  the  concentration  distributions  of  fuel  com¬ 
ponents.  Close  to  the  anode  inlet,  concentrations  of  FI2  and  CO 
increase,  while  concentrations  of  methane  and  FI20  decrease 
resulting  from  steam  reforming  reaction.  As  the  reforming  reaction 
is  enhanced  due  to  H2  consumption  by  electrochemical  reaction, 
methane  is  almost  fully  converted  before  the  middle  section  of 
anode  flow  channel.  At  the  lower  part  of  anode  channel  where 
electrochemical  reactions  dominate,  H2  and  CO  are  reacted  and 
IT20  and  C02  are  produced,  while  the  gas  molar  concentrations 
obey  the  thermodynamic  equilibrium  of  water  gas  shift  reaction. 
Fig.  5(c)  shows  the  polarizations  distribution.  As  fuel  concentration 
decreases  and  cell  temperature  increases,  the  open  circuit  voltage 
appears  to  decrease.  At  the  low-temperature  zone  close  to  the  fuel 
inlet,  the  activation  overpotential  is  relatively  high  and  the  ohmic 
overpotential  possesses  almost  50-60%  in  the  total  overpotential. 
As  the  PEN  temperature  increases,  the  ohmic  overpotential  is  more 
than  70%  of  the  total  overpotential.  So  the  ohmic  loss  is  the  domi¬ 
nant  polarization  in  tubular  SOFCs. 

As  shown  in  Fig.  6(a),  the  local  current  density  shows  a 
non-monotonic  distribution  under  coflow  condition.  The  cell  tem¬ 
perature  is  low  at  the  inlet  section  of  flow  channel,  while  the  fuel 
and  oxygen  concentration  is  low  at  the  outlet  section  of  flow  chan¬ 
nels,  which  leads  to  a  maximum  of  local  current  density  at  the 
middle  section.  Fig.  6(b)  shows  the  fraction  of  FI2  electrochemical 
current  in  the  total  current  density,  where  more  than  80%  current 
density  attributes  by  FI2  electrochemical  oxidation  due  to  higher 
kinetics  of  FI2  than  that  of  CO.  Fig.  6(c)  and  (d)  shows  the  local 
gas  velocities  in  the  anode  and  cathode  channel,  respectively.  With 
a  small  pressure  drop,  the  increment  of  gas  temperature  leads  to 
increasing  gas  velocity  at  the  anode  side,  while  the  gas  velocity  at 
the  cathode  side  shows  a  non-monotonic  variation  due  to  trade-off 
between  the  decrement  of  gas  concentration  and  the  increment  of 
gas  temperature. 

The  neural  network  model  of  turbomachinery  in  Section  2.1.3 
is  validated  separately  based  on  the  performance  maps  of  C/T  from 
Ref.  [37],  which  has  been  used  in  some  literatures  [12,17].  As  shown 


Dimensionless  cell  length 


Fig.  6.  Axial  distribution  of  (a)  local  current  density,  (b)  fraction  of  H2  electrochem¬ 
ical  current  in  the  total  current  density,  (c)  anode  gas  velocity,  and  (d)  cathode  gas 
velocity  in  SOFC. 


in  Fig.  7,  the  neural  network  model  with  preprocessing  of  coor¬ 
dinate  transfer  has  the  ability  to  obtain  a  good  interpolation  and 
extrapolation  accuracy. 


4.  Reverse  design  of  units  and  system  dynamics  analysis 

During  the  part-load  or  dynamic  operation,  performance  of  units 
generally  deviates  from  that  at  design  point,  especially  for  some 
fast-response  units  (e.g.  gas  ejector).  Due  to  lacking  of  open  data, 
reverse  sizing  is  necessary  for  this  rating  problem.  In  this  section, 
based  on  the  full-load  and  stationery  data,  the  recuperator  and  ejec¬ 
tor  are  designed  reversely  for  analysis  of  part-load  and  dynamic 
performance. 


Table  8 

Comparison  of  system  status  among  experiments,  APSAT  and  our  model. 


Status  (T,  p,  rh) 

Experiments 

APSAT 

Our  model 

Compressor  inlet 

15, 1.013,0.635 

15, 1.01,0.635 

15, 1.013,0.635 

Compressor  outlet 

155,3,0.635 

147,3,0.635 

147.2,  3,  0.635 

Recuperator  cool-side  outlet 

500,  2.9,  0.635 

500,  2.94,  0.635 

500.2,  2.9,  0.635 

SOFC  inlet  air  bypass 

500,2.9,  0.14 

500,  2.94,  0.14 

500.2,2.9,  0.14 

SOFC  inlet  air  (after  pipel ) 

500,  2.9,  0.495 

497,  2.94,  0.495 

500.2,  2.9,  0.495 

SOFC  exhaust  (after  burner) 

780,  2.8,  0.502 

770,  2.82,  0.502 

783.06,  2.83,  0.502 

HP  turbine  inlet  (before  pipe2) 

730,  2.8,  0.642 

714,2.82,0.642 

724.7,  2.8,  0.642 

HP  turbine  outlet 

610,  -,  0.642 

575, 1.54,  0.642 

579.8,1.374,  0.642 

LP  turbine  outlet  (before  pipe3) 

550,  -,  0.642 

545, 1.24,  0.642 

551.3, 1.194,  0.642 

System  exhaust 

210, 1.1,0.642 

202, 1.22,0.642 

206.8,1.1,0.642 

SOFC  fuel  inlet 

15,  -,  0.007 

15,2.94,  0.007 

15,2.9,  0.007 
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Fig.  7.  Neural  network  model  of  turbomachinery:  (a)  and  (b)  are  dimensionless  speed  contours  of  compressor,  (c)  and  (d)  are  dimensionless  speed  contours  of  HP  turbine. 
Dashed  lines  are  original  and  solid  lines  are  from  network  outputs. 


4.1.  Reverse  design  of  units 

The  recuperator  here  is  a  cross-flow,  plate-fin  compact  heat 
exchanger,  where  3/8-6.06  blind-type  surface  is  chosen  at  the  air 
side  and  11.11  surface  is  used  at  the  side  of  high-temperature 
exhaust  [29].  Table  9  lists  the  design  parameters  of  recuperator 
and  two  ejectors  from  their  full-load  operating  data. 

As  a  result,  the  recuperator  size  is:  length  of  air-side  flow 
path  is  5.06  m,  length  of  high  temperature  exhaust  flow  path  is 


Table  9 

Design  parameters  of  recuperator  and  ejectors. 


Parameter 

Value 

Recuperator 

Effectiveness 

0.89 

Cold-side  target  relative  pressure  drop  (%) 

0.5 

Hot-side  target  relative  pressure  drop  (%) 

2 

Air  flow  rate  (kg  s-1 ) 

0.635 

Inlet  air  pressure  (bar) 

1.01325 

Inlet  air  temperature  (K) 

419.5 

Exhaust  gas  flow  rate  (kg  s-1 ) 

0.615 

Inlet  exhaust  gas  pressure  (bar) 

1.2 

Inlet  exhaust  gas  temperature  (K) 

817 

Exhaust  composition  (H2O-CO2-O2-N2) 

3.76%  H20-1.89% 

Ejector 

Mass  flow  rate  of  working  fluid  (kg  s-1 ) 

C02-16.85%  02-77.5%  N2 

0.0035 

Working  fluid  temperature  (I<) 

288.15 

Working  fluid  pressure  (bar) 

6 

Mass  flow  rate  of  induced  fluid  (kg  s-1 ) 

0.0186 

Induced  fluid  temperature  (K) 

1300 

Induced  fluid  pressure  (bar) 

2.89 

Induced  fluid  composition 

12.9%  H2-53.5%  H20-10.2% 

Compressed  fluid  pressure  (bar) 

CO-23.2%  C02-0.2%  N2 

2.905 

4.04  m,  and  non-fluid  path  length  is  0.07  m,  which  means  cold- 
side  and  hot-side  frontal  areas  of  0.317  and  0.316  m2,  respectively 
and  volume  of  recuperator  of  1.6  m3.  Note  that  for  a  highly  effi¬ 
cient  gas-gas  compact  heat  exchanger,  due  to  nearly  equal  heat 
capacity  rate  of  hot  or  cold  fluids  in  this  application  (Cmin/Cmax  ~  1), 
the  target  effectiveness  and  target  pressure  drop  has  a  signifi¬ 
cant  impact  on  the  recuperator  size.  For  example,  when  the  target 
effectiveness  is  given  as  £  =  0.85,  the  dimensions  of  recuperator  is 
2.74  x  2.85  x  0.08  m,  which  is  much  smaller  than  that  when  £  =  0.89 
in  Table  9. 

With  the  design  parameter  of  ejectors  in  Table  9,  the  min¬ 
imum  diameter  and  outlet  diameter  of  nozzle  are  2.06  and 
2.07  mm,  respectively,  and  the  diameter  of  mixing  chamber  is 
5.89  cm. 


4.2.  Part-load  and  dynamic  operating 

In  our  multi-level  model  library,  most  of  unit  models  can  be  used 
for  dynamic  analysis.  Fig.  8  shows  a  dynamic  process  corresponding 
to  the  load  variation.  As  shown  in  Fig.  8(a),  due  to  ohmic  loss,  the 
cell  current  decreases  immediately  as  the  cell  voltage  increases, 
then  gradually  goes  to  a  new  steady  state.  The  tendency  of  first 
decrement  and  then  increment  of  the  stack  current  shows  char¬ 
acteristics  of  non-minimum  phase.  Fig.  8(b)  shows  the  variation 
of  fuel  concentrations  at  the  anode  outlet.  As  the  stack  current 
decreases,  the  concentration  of  H2,  CO  and  O2  increase  and  the  con¬ 
centration  of  FI20  and  C02  decrease.  Due  to  the  buffering  effect  of 
volume  of  gas  flow  channel,  the  fuel  concentration  varies  smoothly 
without  a  sudden  change  like  the  cell  voltage-current  relationship. 
Fig.  8(c)  shows  the  transient  injection  coefficient  of  ejectors,  which 
is  defined  as  the  ratio  of  induced  fluid  flow  rate  to  working  fluid 
flow  rate. 
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(a)  Time  /  s 


(°)  Time  /  s 


Fig.  8.  Dynamics  of  (a)  cell  voltage  and  stack  current,  (b)  anode  outlet  fuel  concentration,  (c)  injection  coefficient  of  ejectors,  and  (d)  temporal  variation  and  spatial  distribution 
of  local  current  density  during  a  variable  SOFC  load. 


From  this  part-load  and  dynamic  process,  we  can  find  the  time 
constant  of  SOFC  stack  in  SWPC  demonstration  system  is  almost 
15  s,  which  is  dominant  by  gas  filling  dynamics  in  the  flow  channel. 
As  shown  in  Fig.  6(c),  the  average  velocity  of  anode  gas  is  almost 
0.45 ms-1,  for  the  cell  length  of  1.5 m,  the  corresponding  gas  res¬ 
idence  time  of  3.3  s  is  almost  one  fifth  of  the  stack  time  constant. 
Fig.  8(d)  shows  the  temporal  variation  and  spatial  distribution  of 
local  current  density.  Note  that  the  solid  temperatures  do  not  reach 
the  steady  state  during  the  total  transient  process. 


5.  Conclusion 

This  paper  presented  the  model  libraries  of  turbomachinery, 
steam  methane  reformer,  heat  exchanger,  ejector  and  other  balanc¬ 
ing  units.  The  characteristics  of  our  SOFC-GT  simulation  platform 
are  summarized  as  follows: 

(1)  Compared  to  APSAT  and  other  state-of-the-art  SOFC-GT  simu¬ 
lation  tool,  the  biggest  advantage  is  the  integrated  hierarchical 
modeling  framework  of  natural  gas  IRSOFC,  which  has  been  pre¬ 
sented  in  details  in  Ref.  [2].  And  the  equation-based  solver  of 
gPROMS  is  very  suitable  for  distributed  SOFC  modeling. 


(2)  The  multi-level  models  of  balancing  units  include  part  of  state- 
of-the-art  achievements  and  also  our  specific  work  of,  e.g. 
neural  network  model  with  coordinate  transform  processing 
of  turbomachinery,  uniform  distributed  modeling  of  tube-shell 
and  compact  heat  exchangers,  etc. 

(3)  The  distributed  modeling  is  available  to  reasonable  predic¬ 
tion  of  the  hot  spot  and  the  outlet  gas  temperature  for 
high-temperature  application,  which  is  significant  for  the  safe 
operation  of  SOFC  and  control  strategy  of  TIT. 

(4)  Validated  by  the  steady-state  experimental  data  of  SWPC 
demonstration  SOFC-GT  hybrid  system,  our  prediction  for  sys¬ 
tem  performance  and  system  status  shows  more  accurate  than 
that  from  APSAT  [33]. 

(5)  Most  of  unit  models  can  be  used  for  both  rating  and  sizing  prob¬ 
lem,  and  provide  transient  mechanism  part-load  and  dynamic 
analysis. 

(6)  With  graphic  specification  of  each  unit  model,  our  multi-level 
SOFC-GT  simulation  platform  is  available  for  flexible  configu¬ 
ration  and  modular  design  like  commercial  process  simulation 
software  do. 

All  these  work  provide  a  good  fundament  for  system  design, 

optimization  and  control  algorithms  in  the  future. 
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